% =========================================================================
function outforc=varforc(ydata,lags,nforc,zi,theta) 
% 
% Estimate and forecast using a VAR or BVAR with LAGS 
% 
% Input: 
% -------------
% ydata         [T n]
% lags          integer
% nforc         number of horizons to forecast 
% zi            overall tightness of MN prior
%               set to 0 for OLS 
% theta         rate of decay of MN prior 
%               set to 0 for OLS 
% 
% Output: 
% --------------
% Structure     OUTFORC 
% .By           Autoregressive parameters, Pages correspond to lags 
% .Bx           Constant 
% .u            In sample residuals 
% .yfor         Forecasts [NFORC+1 NZ] 
%               First row is actual data
%
% Notes: 
% 1) Ordering of forecasts 
%    ----------------------
%    First row is actual data 
%
% 2) Parameters of MN prior 
%    -----------------------
% ZI    is overall tightness of the MN prior 
%       Bambura et al use 1/0.26 
% THETA is lag-decay 
%       Bembura et al use 1 
% We abstract here from LAMBDA and MU the own persistece and co persistence
% priors 
% 
% 1. Estimate an OLS VAR with LAGS 
% 2. Forecast NFORC out 
% Alejandro Justiniano May 28 2010 
% =========================================================================
nz=size(ydata,2); 

% -------------------------------------------------------------------------
% 1. VAR 
% -------------------------------------------------------------------------
if (zi==0 && theta==0)
    disp('Classical OLS VAR') 
    outforc=rfvar3(ydata,lags,ones(length(ydata),1),[],0,0);
else 
    disp('Bayesian VAR'); 
    [By,Bx,u,xxi,y,X]=rfvar3_min_forLSTVVAR(ydata,lags,ones(length(ydata),1),[],[],[],zi,theta); 
    outforc.By=By;  
    outforc.Bx=Bx;  
    outforc.u =u; 
end 
outforc.By=reshape( outforc.By,[nz nz*lags] ); 

% -------------------------------------------------------------------------
% 2. Forecast 
% Obtain companion representation of VAR(p) for p > 1 
% Forecast NFORC+1 periods ahead 
% -------------------------------------------------------------------------
if lags > 1
    PHI=zeros(nz*lags,nz*lags);
    PHI(1:nz,:)=reshape( outforc.By ,[nz nz*lags] );
    PHI(nz+1:end,1:nz*(lags-1) ) = eye(nz*(lags-1) );
    PHIC=[outforc.Bx(:);zeros(nz*(lags-1),1)];
else
    PHI=outforc.By;
    PHIC=outforc.Bx(:);
end

yfor     =zeros(nforc+1,nz*lags ); 
yfor(1,:)=reshape( flipud( ydata(end-lags+1:end,:) )' , [nz*lags 1] ); 
for ii=2:nforc+1; 
    yfor(ii,:)=(PHI*yfor(ii-1,:)'+PHIC); 
end 
outforc.yforc=yfor(:,1:nz); 